Particle-Number Restoration within the Energy Density Functional formalism: 
Nonviability of terms depending on noninteger powers of the density matrices 

T. Duguet,i'2,3,g ]y[ Bender,-*. 5.0 K. Bennaceur,6.3.[i| D. LacroixJ.I and T. Lesinski^.^ 

National Superconducting Cyclotron Laboratory, 1 Cyclotron Laboratory, East-Lansing, MI 48824, USA 
^Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA 
'^CEA, Centre de Saclay, IRFU/Service de Physique Nuclaire, F-91191 Gif-sur-Yvette, France 
^ Universite Bordeaux, Centre d'Etudes Nucleaires de Bordeaux Gradignan, UMR5797, F-33175 Gradignan, France 
^CNRS/IN2P3, Centre d'Etudes Nucleaires de Bordeaux Gradignan, UMR5797, F-33175 Gradignan, France 
^ Universite de Lyon, F-69003 Lyon, France; Universite Lyon 1, F-69622 Villeurbanne, France; 
CNRS/IN2P3, UMR 5822, Institut de Physique Nucleaire de Lyon 
''GANIL, CEA et IN2P3, BP 5027, 14076 Caen Cedex, France 
(Dated: March 4, 2009) 

We discuss the origin of pathological behaviors that have been recently identified in particle- 
number-restoration calculations performed within the nuclear energy density functional framework. 
A regularization method that removes the problematic terms from the multi-reference energy den- 
sity functional and which applies (i) to any symmetry restoration- and/or generator-coordinate- 
method-based configuration mixing calculation and (ii) to energy density functionals depending 
only on integer powers of the density matrices, was proposed in [D. Lacroix, T. Duguet, M. Ben- 
der, arXiv: 0809. 2041 and implemented for particle-number restoration calculations in [M. Bender, 
T. Duguet, D. Lacroix, arXiv: 0809. 2045 . In the present paper, we address the viability of non- 
integer powers of the density matrices in the nuclear energy density functional. Our discussion 
builds upon the analysis already carried out in [J. Dobaczewski et ai, Phys. Rev. C 76, 054315 
(2007)]. First, we propose to reduce the pathological nature of terms depending on a non-integer 
power of the density matrices by regularizing the fraction that relates to the integer part of the expo- 
nent using the method proposed in [D. Lacroix, T. Duguet, M. Bender, arXiv:0809.2041 . Then, we 
discuss the spurious features brought about by the remaining fractional power. Finally, we conclude 
that non-integer powers of the density matrices are not viable and should be avoided in the first 
place when constructing nuclear energy density functionals that are eventually meant to be used in 
multi-reference calculations. 

PACS numbers: 21.10.Re, 21.60.Ev, 71.15.Mb 
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I. INTRODUCTION 

In their recent paper Dobaczewski et al. have 
pointed out that there are two distinct pathologies that 
might appear in calculations aiming at restoring particle 
number within the nuclear Energy Density Functional 
(EDF) framework. Formulating a Particle Number Re- 
stored (PNR) EDF calculation through a contour integral 
in the complex plane over multi-reference (MR) EDF ker- 
nels, the two categories of pathologies are associated with 
spurious poles and branch cuts of the complex MR-EDF 
kernels that relate to dependencies of the latter on in- 
teger and non-integer powers of the (transition) density 
matrices, respectively. 

The possible appearance of spurious poles was already 
identified in Refs. 0, [1, Q. In Ref. hereafter re- 
ferred to as Paper I, we demonstrated that such a pathol- 
ogy is shared by any symmetry-restoration- or generator- 
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coordinate-method (GCM)-based configuration mixing 
calculation performed within the EDF context, which 
we will call a Multi- Reference Energy Density Functional 
(MR-EDF) formalism from hereon. In most other cases 
than PNR, however, the identification of the spuriosities 
is much less transparent. In Paper I, we proposed a for- 
mal and practical regularization method that applies to 
any symmetry restoration and/or GCM-based configura- 
tion mixing calculation. In Ref. hereafter referred to 
as Paper II, we applied the correction method to PNR 
calculations using a particular energy functional that de- 
pends only on integer powers of the density matrices and 
thus only display spurious poles. 

The pathology associated with spurious branch cuts 
has been overlooked until recently [l| for reasons that 
will become clear in the following. As a remedy to it, the 
authors of Ref. [l| have proposed to deform the integra- 
tion contour in the complex plane such that it does not 
cross the branch cuts. As will be discussed below, such 
a procedure does not allow the definition of a fully sat- 
isfactory theory; e.g. the breaking of the shift invariance 
remains. In addition, there is no clear method for gen- 
eralizing the proposed solution to any other coordinate 
frequently used in MR-EDF calculations. 

In the present paper, we thus address the pathology as- 
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sociated with branch cuts from a different point of view 
than in Ref. We first make use of the correction 
scheme designed in Paper I to regularize the pathology 
associated with spurious poles. Doing so we can isolate 
the part that is specific to the pathology brought about 
by branch cuts and question whether it is possible to 
perform meaningful MR calculations using an EDF that 
depends on non- integer powers of the density matrices. 
In fact, the question relates to the possibility to deal 
with any EDF providing multi-valued MR kernels over 
the complex plane. It will appear that any EDF (i) pro- 
viding multi- valued MR kernels over the complex plane, 
(ii) whose functional form is such that the pole structure 
cannot be extracted analytically; e.g. the family of func- 
tional proposed by Fayans and collaborators 0, @1, is 
critical. Eventually, anything but low-order polynomials 
seems difficult, if not impossible, to handle in practical 
MR-EDF calculations. Indeed, even if the pole structure 
of a complicated EDF can be characterized, it is only 
for low-order polynomials that the regularization method 
proposed in Paper I can be applied to identify the associ- 
ated spurious contribution to the physical pole at z = 0. 

The present discussion is conducted for PNR calcula- 
tions based on a EDF whose normal part takes the form 
of a toy Skyrme energy density functional, and whose 
pairing part derives from a density-dependent delta in- 
teraction (DDDI). Numerical applications arc performed 
using the realistic SLy4 Skyrme EDF combined with a lo- 
cal pairing part as derived from a (density-independent) 
delta interaction (DI). Two situations of interest are ac- 
tually considered that correspond to using an EDF (i) 
derived from (density-dependent) forces (ii) formulated 
directly at the level of the energy functional itself. 

The paper is organized as follows. In Scct.|ITl basic el- 
ements of the single-reference EDF method are recalled 
and the form of the simplified energy functional consid- 
ered for the discussion is given. Section IIIII introduces 
PNR calculations performed within the EDF framework 
and describes the analytical continuation into the com- 
plex plane which is used for analysis purposes in Sect. lIVl 

Section lTVl discusses the occurrence of pathological pat- 
terns in particle-number restored energies. First, we re- 
call the situation for EDFs depending on integer pow- 
ers of the density matrices, which is the focus of Pa- 
pers I and II. Then, EDFs depending on non-integer pow- 
ers of the density matrices are discussed as the simplest 
and most practically relevant example of EDF generating 
multivalued PNR energy kernels over the complex plane. 
Still, the conclusions drawn are valid for more involved 
EDFs presenting such a feature. Finally, results of nu- 
merical applications are provided in Sect. El highlighting 
again the differences between EDFs depending on integer 
powers of the density matrices and those depending on 
non-integer ones. Conclusions are given in Sect. IVIl 



II. SINGLE-REFERENCE EDF METHOD 

Before we present results obtained with a realistic 
SLy4-|-DI EDF, we analyze the relevant physics with a 
toy functional, reduced to the bare minimum of terms 
necessary to convey our point. 

A. Density Matrices 

The implementation of the Single-Reference EDF ap- 
proach relies on the use of a quasi-particlc vacuum \^;p) 
to calculate the one-body density matrices the energy 
£[p,K,K*] is a functional of. The index in |3><p) de- 
notes the gauge angle that provides the orientation of 
the system in gauge space. Using the requirement that a 
meaningful energy functional should be invariant under 
gauge space rotations, the angle can be set to a conve- 
nient value, usually (f = 0. 

In the canonical basis {(/>^(r) = (rjaJ^lO)} of the Bo- 
goliubov transformation that underlies the quasi-particle 
vacuum |$o), the SR normal density matrix p and 
anomalous density matrix k (pairing tensor) take the 
form 



($o|ata^|^n) 

(*o|*o) 
($o|a^a^|i>o) 

(*o|*o) 
(^ol4"tl^o) 
($o|«'o) 



^1 Sf,^ , 



(1) 
(2) 

(3) 



where are BCS-like occupation numbers such 



that uj, + vz 



1, = M 



^ > and 



-Vfi. The two 



canonical states (/i, p) are the so-called pair conjugated 
states. Based on an appropriate quantum number, the 
basis can be split into a positive half {p > 0) and a 
negative half {p < 0). When a canonical state iJ. belongs 
to one of these halves, its conjugate state jl belongs to 
the other half. 

From the point of view of their physical content, cur- 
rently used nuclear EDFs can be put under the generic 
form 



£[p, K, K*] = fkin[p] + ^norm[p] + £pa.iT[p, «, K* 



(4) 



where appear the uncorrelated kinetic energy, the normal 
and the pairing contributions, respectively. The contri- 
butions from Coulomb interaction and explicit quantum 
corrections as the center of mass correction have been 
omitted for the sake of using simple notations. Including 
them will not modify the arguments given below. From 
the point of view of their functional dependence on nor- 
mal and anomalous density matrices, the different parts 
of the functional can be formally written 



(5) 
(6) 
(7) 
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where the superscripts specify the powers of the nor- 
mal and anomalous density matrices that contribute to 
a given term. The focus of the present work is on the 
properties of Sppp" and when such a nuclear EDF 

is used in MR calculations. Note that the Slater approxi- 
mation which is usually used to tackle the exchange part 
of the Coulomb contribution to the normal part of the 
EDF is of the form Eppp"^ . 

For the sake of a transparent discussion, we will per- 
form the analysis for a toy functional limited to the min- 
imum of ingredients necessary to make the point. For 
this purpose, we start from a simplified Skyrme interac- 
tion containing the so-called and ^3 terms only, which 
limits the local densities entering fnormip] to those that 
do not contain spatial derivatives 0]. In addition, and 
because the validity of the points made together with the 
conclusions reached do not critically depend on it, we 
omit the isospin degree of freedom and consider one nu- 
cleon species only throughout the discussion. Comments 
on the additional complexity brought by considering neu- 
trons and protons are added in Sect. IIVEI The gener- 
alization to a complete and realistic Skyrme or Gogny 
EDF is then straightforward. 

B. Local densities 

The local matter and spin densities needed to construct 
£nornAp] ^re givcu by 

p{v) ^ ^</.t(r)0^(r)p^^, (8) 

s(r) = ^0j,(r)o-0^(r)p^^, (9) 

where 4>^{r) and & denote a canonical single-particle 
spinor and the vector of Pauli matrices, respectively. In 
addition, one needs the local kinetic density 

T(r) ^ 5^[V0t(r)] . [V0^(r)]p^^, (10) 

to express the kinetic energy. The three previous local 
densities can be put under the form 

/(r) = ^Ty^(r)p^^, (11) 

where / e {p, s, r} and where the explicit form of Wj^^{r) 
can be easily extracted from Eqs. (|5]fTU)): i.e. 

M/;,(r) EE 0t(r)0.(r), (12) 

W^,(r) EE 0t(r)^0,(r), (13) 

W;,{v) EE [V0t (r)] . [V0.(r)] . (14) 

The densities entering the pairing part of the EDF arc 
the local pair densities defined as 

p(r) EE 2^<^(r)«^^. (15) 

M>0 



Finally, with the symmetries of the SR and MR EDF 
calculations assumed here, W^p(r) and W^^ are equal 
and given by the spin-singlet part of the two-body wave 
function, defined as 

W^;.(r) = M^;;(r) ^ J2 '^Mr^)Mr-^) (16) 
(j=±i 

= -WP^{r) = ~W^;{v). (17) 

C. Toy energy density functional 

The kinetic energy part of the EDF takes the form 

8^ EE /rfV|lr(r), (18) 

whereas the normal part derives from a toy Skyrme in- 
teraction characterized by^ 

£PP = J d^r[APPp'^{r)+ A"' s'^{r)] , (19) 

SPPP" EE jSr [APPP" p2 (r) + A^'P" (r)] p" (r) .(20) 

Finally, the pairing part of the EDF is given as 

^ ld\A^Pp*{r)p{r), (21) 

e'^^'P'' EE j cPrAP~PP' p*{r)p{Y)p^{r), (22) 

where the superscripts // and ///' of the As refer to the 
local densities the corresponding term depends on. In ad- 
dition, one can still read off those superscripts the powers 
of normal and anomalous density matrices that the cor- 
responding term incorporate. Note that no hypothesis 
about time-reversal invariance of the system has been 
made. On the other hand, we limit ourselves to quasi- 
particle vacua 1$^} with an even number-parity quan- 
tum number and thus only discuss explicitly even-even 
systems. 

The part of the EDF which only depends on the normal 
density matrix can be derived from a schematic Skyrme 
force 

?;sfc(R, ri2) = to{\ + xaPa)5{'^i2) 

+ |(l + x3P,)p[?(R)<5(ri2), (23) 

where R = (ri + r2)/2 and Y12 = ri — r2, whereas 
Pa = \{1 + cTi ■ (J2) denotes the spin exchange operator. 
Computing the normal part of the EDF as the Hartree 



^ One could have considered that the terms multiplying and s'^ 
in £PPP present different exponents. 
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and Fock contributions derived from such an empirical 
effective vertex, one obtains 



= +lto{l - x„) , APPP" 



-l^t^{l-x^), (24a) 
-^t3(l-X3), (24b) 



which shows that in this case the four couphng constants 
entering the EDF depend on two independent parame- 
ters only. However, we will also be interested in EDFs 
which are not derived from a Skyrme force and for which 
the four coupling constants can be chosen independently. 
For more complete and realistic functionals, local gauge 
invariance imposes constraints between certain coupling 
constants [loj . 

The part of the EDF which depends on the anomalous 
density matrix could be derived from the same Skyrme 
force. As one usually focuses on the superfluidity in the 
spin-singlet /isospin-triplet channel, one would be led in 
practice to select only a part of the interaction from the 
outset. Furthermore, there exists strong theoretical mo- 
tivations to explicitly disconnect the part of the EDF 
responsible for superfluidity from t he p art that only de- 
pends on the normal density matrix . However, such a 
decoupling between f^norm and £pair is at the origin of seri- 
ous problems encountered in MR-EDF calculations P, Q . 
We will come back to that in the following. For now, one 
can relate the specific local pairing functional given in 
Eqs. (|21j|22p to a DDDI vertex of the form 



i;pair(R,r)-|(l-P.) 



Po(R)\^ 



Psat 



<5(r), (25) 



where psat — 0.16 fm which leads to 



APP = \i,., 



^io. (26) 



Independently of the starting point, a quasi-local pair- 
ing EDF must be regularized/renormalized as far as its 
ultraviolet divergence is concerned (l^ . 



III. PARTICLE NUMBER RESTORATION 
A. Notations 



in such a way that £^ depends only implicitly on the 
(normalized) projected state 



($o|P^|$o) 



2tt 



dip 



-iipN 



2tt cn 



\%)- (29) 



The gauge-space-rotated product states constituting the 
MR set of interest read, in their common canonical basis, 



as 



1$. 



n 

p>0 



2iip + + 



|0), (30) 



where |0) is the particle vacuum. The above form of 

is convenient to compute the overlap between a rotated 

state and the unrotated one 



Hi 

li>0 



(31) 



In Eq. ([271), 5[0, denotes the (set of) MR energy den- 
sity functional kernel(s). It is traditionally defined by 
replacing the SR normal and anomalous density matri- 
ces by transition ones 















(*o|<f>^) 




($o|ai4l*v> 




('&o|<f>^) 





(32) 
(33) 



5.-^. (34) 



into the SR EDF £[p, k, k*]. This corresponds to defining 
non-diagonal energy kernels through the prescription 



(35) 



As discussed in Paper I, MR-EDF calculations per- 
formed along the lines presented above fulfill basic re- 
quirements [13| but may display pathologies such as di- 
vergences and finite steps in the energy. The extent of 
such problems depends on the analytical form of the EDF 
used. In order to conduct an in-depth analysis of the po- 
tential problems, it is necessary to perform an analytical 
continuation of £[0, ip] to the complex plane [ll. [T^. 



Continuation to the complex plane 



As extensively discussed in Ref. [l| and in Paper II, 
Particle Number Restoration (PNR) performed within 
the EDF framework relies on calculating the energy of 
the A'^-particle system through a MR energy functional 
of the form 



The continuation is achieved by extending the complex 
number z = e^'P onto the entire complex plane in all 
previous formulae.^ In that context, the PNR energy 



where 



n2TT —iipN 

/ dip—^£[0,p](<i>o\%), (27) 

Jo ^TTC^ 



c% = 



2ir ^-iNip 



dp ■ 



2-K 



(28) 



^ The same notation as before is used when extending the defini- 
tion of SR states and energy kernels to any value of the complex 
variable z. Thus, we abusively replace the gauge angle ip by the 
complex variable z in all our expressions; i.e. SR states character- 
ized by the gauge angle i^, |3>^) are extended into l^^) to denote 
SR states anywhere on the complex plane. In particular, the un- 
rotated SR state, denoted as |<E>o) when using as a variable, is 
written as when using z as a more general variable. 
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defined through Eq. p7|) results from integrating over 
over a closed contour around z = which can be chosen 
as the unit circle Ci (|z| = R = 1) 



dz £ [z\ 



C, 2iTTC% zN+ 



1 



dz 



2iTT 



where 



M>0 



(36) 
(37) 

(38) 



With this continuation, the transition density matrix and 
pairing tensor read as 





9 2 








K^"" - 


Uf^ Vf^ z 


2 


ul + -"l 


z2 


zl * 






ul + vl 


Z2 



(39) 

(40) 
(41) 



and must replace the SR density matrices in Eas.[51 fT51 in 
order to define the corresponding transition local densi- 
ties. Finally, the energy kernel from Eq. ([55]) reads as 



E[z]=£[p^\k^\k'-^*] 



(42) 



For PNR calculations, this simply amounts to expressing 
the EDF kernel £[z\ in the canonical basis of the Bogoli- 
ubov transformation defining any of the product states of 

t 

Im[z] V +Z3 



Q +Z2 

I 
I 

I 

9 +Zi 

I 

— I 

z, 



Z= i lu^l/lv^l 



o 







Re[z] 



FIG. 1: Pole structure of 
plane. 



and * on the complex 



reference; e.g. Indeed, the same canonical basis is 

shared by all product states 1$^) over the complex plane, 
as well as by the Bogoliubov transformation linking any 
pair of them. 



IV. STEPS AND DIVERGENCES 

A. General considerations 

The computation of £^ through an integration over 
a contour encircling the origin requires the knowl- 
edge of the (non-)analytical structure of the integrand 
£[z] ($i|<i>2)/z^+^ over the complex plane. First, it ob- 
viously contains a (physical) pole at z = 0. Since £[z\ 
is a functional of the transition density matrices, (i) it is 
a function of and is thus even, i.e. £[z\ = f [— z], (ii) 
its analytical structure relates to the one of the tran- 
sition densities. As displayed in Fig. [U it is trivial 
to see that p^^ , k^^ and k^^* possess simple poles at 
z = ±z^ = ±i|up|/|w^| [11]. In general, it is likely that 
those poles will translate into non-analytical features of 
£[z] ($i|<i>z} that have serious consequences on the PNR 
energy. 

As explained in Paper I, it is necessary to go to con- 
figuration space to isolate the spurious contributions to 
the MR-EDF energy. For a given pair of vacua belong- 
ing to the MR set, the basis relevant to the analysis of 
the corresponding energy kernel is the canonical basis of 
the Bogoliubov transformation connecting the two vacua. 



B. Term depending on integer powers of densities 

Let us start the analysis with terms that depend on 
integer powers of the density matrices. To illustrate the 
situation, we make use of the bilinear parts, Eqs. (jl9p 
and dH]), of the toy EDF introduced in Sect. HTC] 



1. Matrix elements 

Working in the canonical basis of Bogoliubov transfor- 
mation connecting |$i) and the bilinear part of the 
energy kernel £[z] takes the form 

£PP[z] + £''^[z] 

where v^^ and v'^'^ denote matrix elements of effective 
two-body vertices associated with £pp and £'^'^, respec- 
tively. For the toy functional of Eq. (|19p . the matrix 
elements of v^p take the form 
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^2jd\ [APP WP^{v) WP^ir) + W^^(r) • W:,(r)] 



(44) 



r 



The quasi-local nature of the Skyrme energy functional 
(the toy functional considered here being purely local) 
simplifies the construction of the matrix elements v^f^^^, 
as they involve a single spatial integral only. However, 
the discussion conducted in the rest of the paper would 
hold equally for non-local functionals; e.g. as obtained 
from finite-range, possibly non-local, effective vertices. 

The matrix elements associated with f in Eq. ((2T|) 
take the form 



ing on integer powers of the density matrices do not de- 
pend on the pair of vacua |<I>i) and under consider- 
ation, i.e. they do not depend on the gauge variable z. 



2. Analytical structure of {S''''[z] + £'^'^[z\) {$i|$^> 



^Ajd\ APP WP; (r) WP, (r) . (45) 

Note that for PNR calculations, the matrix elements that 
one naturally associate to any term of the EDF depend- 



Due to the additional presence of the norm factor 
('i'll'i'z) in the integrand of Eq. ([55]) . it is easy to realize 
that only the terms corresponding to v — ^ and u = fi in 
Eq. can lead to non-analytical features 0, Q . Such 
terms contribute to the integrand through 



J 



^(vP 

2 



fP 
2 \ MMMM 



4 4 



.22 



n 



fj,fj,IJ,fi fi 



2 9 



v^z^ 



n 



(46) 
(47) 



r 



and both contain potential poles at 2 = = 
±i|u^|/|t;^|. Note that those poles do not exist in the 
first place if the states (^, p.) arc more than doubly de- 
generate in terms of occupation numbers as an additional 
factor from the norm then compensates the single pole 
in Eqs. (gUTlj.^ 

Otherwise, the poles disappear in Eq. (jiGj) if, and only 
if, yPP 



-nn 

7;- - - - 



0; i.e. the matrix elements associ- 
ated with £PP are antisymmetrizcd. Coming back to the 
toy Skyrme functional used in the present paper, and 
noticing that 



|W^,(r)p = |W^;^(r)p= [ ^ l^^(rcT) 

a=±l 



(48) 



for all /i, one finds that vPP 
if. A"" = 



-nn 

V-- - - 



MMMM - -MMMM - if, and only 
APP. As shown by Eqs. (|24all24bp . such a 



^ This holds for bilinear functionals. A term of order n in the den- 
sity matrices can generate a pole at it^^j of order (at most) (n— 1). 
For the pole to disappear, (n— 1) additional factors from the norm 
kernel are needed to cancel the denominator (n^ -l-''^ ^2^ — (n— 1) 
Thus, the pair of interest {p., p.) needs to be degenerate (at least) 
with (n — 1) other pairs in terms of occupations for this to occur. 



condition is satisfied when starting from the (density- 
independent part of the) Skyrme force. The previous 
analysis is trivially extended to the density-independent 
part of a more complete Skyrme or Gogny vertex. On the 
other hand, using a functional approach that bypasses 
the introduction of a two-body vertex, relationships such 
as A'"' = —APP might not be fulfilled. In such a case £pp 
generates poles at z = in the integrand of Eq. ((36)) . 

The poles disappear from Eq. (|T7|) if, and only if, 
^%.pLfi = ^'"ilfLup,! diagonal matrix elements involving 
two conjugated canonical states are identical in £pp and 
f*". If it is so, the two terms in the bracket of Eq. (jTZ]) 
combine in such a way that the dangerous denomina- 
tor explicitly cancels out. One is then left with a finite 
contribution to the MR energy kernel. Such a recom- 
bination is obviously satisfied if both £pp and f are 
constructed from the same (effective) force, for example 
when using the density-independent part of the Gogny 
interaction Using a functional approach or starting 
from two different effective vertices to build £pp and f 
the recombination is unlikely to occur and one is left with 
an ill-defined PNR formalism and compromised results. 



Just as we did to ensure that 



r,PP 



0, i.e. 



A^'^ = —APP, one could work out minimal constraints 
between the coupling constants entering £pp and f "'^ to 
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Im[z] 



c =c +c 





/ ] 












\ \ 




Re[z] 



small radii. Using such contours, it is easy to prove that 



<?N 



Ties 
TZes 



z=0 



z 



N+1 



(49) 
(50) 



2=0 



Because the only pole of the integrand is at z = 0, the 
same result is obtained for £^ by starting from any inte- 
gration contour encircling the origin in Eq. (j36p . When 
the energy is calculated as the average value of a Hamil- 
tonian in the projected state, the independence of the 
projected energy on the details of the integration con- 
tour, as for example its radius, can be related to the 
invariance of the normalized projected state with respect 
to shift transformations P, [l^ . This symmetry will be 
discussed below in the EDF context. 



FIG. 2: Computation of £^ for an EDF (i) obtained from 
the average value of a genuine Hamiltonian in the projected 
state (ii) depending only on integer powers of the densities 
and after applying the correction proposed in Paper I. The 
integration is performed in the complex plane over a circular 
contour Cr of arbitrary radius R. 



4- PNR energy from an EDF 



impose that wf^'^^n 



in the underlying EDF. 



3. Projected energy from a Hamiltonian 

As seen from the previous discussion, poles in the 
transition densities do not always translate into poles 
in £[z] {^i\^z)- The most trivial example for this oc- 
curs when the particle number projected energy is com- 
puted from the average value of a genuine Hamiltonian 
in the projected state l^*^}; i.e. what we denote as the 
strict projected HFB approach in Paper II. In this case, 
the only pole of the integrand in Eq. (j36|) is the phys- 
ical one at z = 0. To apply the Cauchy theorem"' 
and calculate the projected energy, the original circu- 
lar contour Ci must be deformed to exclude the pole at 
z = 0. As shown in Fig. [2l this can be achieved by 
choosing two semi-circular contours Cd and Cg, such that 
Ci = [Cg-I-Crf] (e — > 0), and by closing those semi-circular 
contours along the imaginary axis in such a way that the 
pole at z = is bypassed by two semi-circles of infinitely 



The present Section reformulates parts of the analysis proposed 
in Paper II for functionals proportional to integer powers of the 
density matrices, i.e. we employ Cauehy's integral theorem rather 
than using directly Cauehy's residue formula. Coming back to 
Cauehy's integral theorem will be needed to conduct the dis- 
cussion for more general functionals as is indicated in the next 
Section. 



The poles subsist in Eqs. 



and (gZl) for any EDF 



that is characterized by vfP,^^ ^ and/or ^ }'^'^f^p.- 

To apply the Cauchy theorem in this case, the circular 
contour Ci must now be deformed to exclude not only 
the pole at z = but also those at z = ±z^ which are 
inside the unit circle. As shown in Fig. [3l this can be 
done by choosing two semi-circular contours Cd and Cg , 
such that Ci = [Cg + Cd] (e — * 0), and by closing each 
of them along the imaginary axis in such a way that all 
the poles are bypassed by semi-circles of infinitely small 
radii. Using such contours, the Cauchy theorem leads to 



N 



E 

=0,±z„ 



TZes 



(51) 



whereas remains unchanged. 

According to Eq. ([5T|) . the existence of poles at z = 
±z^ in the integrand makes the PNR energy to (i) depend 
on the radius of the integration circle P, Q (ii) display a 
finite step whenever a pole leaves the integration circle; 
e.g. as the system is deformed along a collective degree 
of freedom Q . Such a behavior make the PNR energy 
to break shift invariance. This is very undesirable as the 
concept of shift transformation and shift invariance can 
be extended to the EDF framework in such a way that 
the invariance of £^ with respect to the radius of the 
integration contour remains a fundamental feature of the 
theory [H. 

Also, PNR energies may display divergences whenever 
a pole crosses the integration circle. When a pole sits on 
the integration contour Cn, the definition of the contour 
Cn. = \Cg + Cd\ (e — *■ 0) is in fact ambiguous and requires 
an additional prescription. The most natural procedure 
is to define the integration through the pole in the sense 
of the Cauchy principal value. Doing so provides a finite 
PNR energy if the Laurent series of the integrand cen- 
tered at the pole only contains odd powers. Considering 
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im[z] I c,=q+c. 




Re[z] 



FIG. 3: Computation of £ for an EDF depending on integer 
powers of the densities. The integration is performed in the 
complex plane over the unit circle Ci. 



the structure of the nuclear EDF. this will happen if the 
EDF (i) only contains bilinear terms (ii) contains addi- 
tional trilinear terms that do not allow three powers of 
the same isospin (as a zero-range three-body force does 
not allow) (iii) contains additional quartic terms which 
are bilinear in each isospin. In this case, one is left with 
simple poles at z — ±z,, and the Cauchy principal value 
equals half the result that would be obtained if the pole 
were to lie inside the integration circle. In all other cases, 
one can see that (i) the poles at z = ±z^ will be of higher 
orders (ii) the Laurent series centered at those poles will 
contain even powers (ii) the Cauchy principle value will 
lead to an infinite values and the PNR energy will di- 
verge as a poles crosses the integration circle. If the EDF 
used is such that PNR energies diverge whenever a pole 
crosses the integration circle, it is important to note that 
Variation After Projection (VAP) calculations will not 
converge as soon as the minimization procedure "finds" 
the infinity 

All previous features prove that PNR calculations are 
ill-defined whenever poles at z ^ arise and that the 
theory is unacceptable as it is. However, it is possible to 
meaningfully regularize PNR calculations based on any 
EDF depending on integer powers of the density matri- 
ces as was demonstrated in Paper I and exemplified in 
Paper II. As a matter of fact, the method proposed in 
Paper I precisely removes the poles at z = ±z^ from 
£[z] (<I>i|$z). However, it is crucial to realize that the 
correction method does not only remove those poles but 
also consistently subtracts a spurious contribution to the 
physical pole at z = 0]. In the end, only the physical 
pole at z = remains in Eq. P5|) and the independence 
of £^ on the integration contour is recovered, as seen 



from Fig. [2l i.e. the same PNR energy is obtained by 
integrating over circular contours Cr of arbitrary radius 
R. 



C. Non-integer power of densities 

1. Problem 

The situation is often more complex due to the pres- 
ence of higher-order terms of the form £ppp and £''^'^p 
in realistic nuclear EDFs, Eqs. and (|22)) . 

If a = 7 = 1, then £ppp and £'^''^p can, at least formally, 
be analyzed as if they originated from a three-body ver- 
tex. Thus, and as for the bilinear terms, two cases have 
to be distinguished (i) and are both derived 

from the same antisymmetrized three-body vertex and 
do not lead to divergences and steps in MR-EDF calcu- 
lations (ii) they refer to different three-body vertices such 
that the regularization method proposed in Paper I can 
be applied to obtain a meaningful PNR-EDF method. 

However, all modern parameterizations of the nuclear 
EDF, starting either from a functional approach or from a 
density-dependent vertex, depend on non-integer powers 
of the density matrix that one cannot expand in a Taylor 
series to relate them, at least formally, to three-body, 
four-body, . . . forces. The goal of the present paper is 
to characterize the pathologies brought about by such 
dependencies and whether or not they are viable in the 
end; i.e. if the corresponding pathologies can be easily 
regularized. 

2. Regularizing the integer part 

As a first step, one can reduce the extent of the prob- 
lems associated with terms of the form £'p^'"+"+'' ^nd 
£^ P , with m and n integer, and < a < 1 and 
< 7 < 1, to pathologies only due to the fractional 
powers and , respectively. This means that steps 
and potential divergences associated with the integer part 
2m + n can be regularized from the outset. This is the 
case either (i) if one started from a density-dependent 
{2m + n)-body effective force or (ii) by applying the 
correction method proposed in Paper 1 to £p ^ and 

Let us exemplify how an empirical extension of the 
correction method proposed in Paper I can be designed 
to regularize the quadratic part of £ppp , with < a < 1. 
To simplify the situation further, we disregard the term 
gKKp ^j^g following discussion. Such a simplification 
does not alter any of the conclusions given in the rest of 
the paper. 

To proceed, we first introduce pseudo two-body ma- 
trix elements vPPf^j^[z] which take, for the toy functional 
considered in the present paper, the form 
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[z]^2jd\ [A"""" M^;^ (r) (r) + A^^"'' W^^ (r) ■ W^, (r)] [p'^ (r)] " . (52) 



With the pseudo two-body matrix elements w^C^^I^] 
at hand, one can apply the correction formula given by 
Eq. (43) of Paper II. However, and as opposed to terms 
of the EDF depending on integer powers of the density 
matrices, the matrix elements of v'^i'f do depend on the 
gauge variable z. As a result, Eq. (43) of Paper II must 
be applied in such a way that the matrix elements are 
located underneath the integral over z. Last but not 
least, it would also be trivial to regularize the integer 
part of E'^'^P by introducing the pseudo two-body ma- 
trix elements v'^^f [z\ and by using them in Eq. (43) of 
Paper II. 

D. Left-over fractional power 

With the latter correction at hand, the quadratic part 
of does not create any divergence or step in the 



where 



with N even. In agreement with the properties of £[z\ 
mentioned above, F[z](r) is an even function of z for all 
r. For odd A^, it is easy to prove that F[z](r) is an odd 
function of z in such a way that F\z\{y) / z-^ remains 
itself an odd function of z. 

The terms corresponding to v = p and v = p are 
absent in Eq. ([M]) because (i) they were removed by the 
correction method briefly outhned in Sect. IIVCI (ii) one 
started from a density-dependent two-body interaction; 
i.e. the term with i' — p do disappear {A^'^p = —Appp ) 
whereas the term with u ~ p could be combined with 
the corresponding one in E'^'^i^ to give a well-behaved 
contribution that we omit here. 

To understand the features displayed by the contri- 
bution S'^Ippp"] to the PNR energy, it is necessary to 
extract for each r the non-analytical structure of the in- 
tegrand in Eq. ([55]) where the order of the two integrals 
over r and z have been reversed. Clearly, the function 
F[z]{y) / z'^ displays a (physical) pole at z = 0. The 



I 

PNR-EDF energy anymore. Again, the same is true if 
one starts from the outset from a density-dependent two- 
body antisymmetrized interaction, as long as the corre- 
sponding term S'^'^'' is explicitly considered in the EDF 
to proceed to the necessary recombination of terms in 
Eq. ((47|) . One way or another, one is only left in the end 
with discussing the impact of the fractional power of the 
transition density; i.e. the extra factor [p^^(r)] , with 
< a < 1. 



1. Analytical structure of £'''''' [z] (<&i |<1?2) 

Now that the pathologies due to the bilinear factor in 
£PPP have been taken care of, the contribution of interest 
to the PNR energy can be written as 



(53) 



(54) 

I 

difficulty comes from the fractional power of the local 
transition density that multiplies i^[2](r). Indeed, such a 
function is multivalued on the complex plane for all r. 

Defining the function corresponding to taking the frac- 
tional power of a complex number^ requires the introduc- 
tion of a branch cut along the axis where that number 
is real and negative. Here, this means that one needs 
the values of z for which the function p^^(r) is real and 
negative. As can be seen from Eqs. ([5]) and ([M)) . the 
transition density is real both on the real and imagi- 
nary axis, but can be negative only on the latter. A 
discussed in Ref. p^^(r) is negative for z = iy such 



^ Parameterizing z = re'*, 9 € [— vr, +7r], we define the principal 
value of the function z" , a being a rational number between zero 
and one, as z"' = r°e"^*. The latter choice lifts the ambiguity 
regarding the multivalued nature of the function but requires to 
track the latter through several Riemann cuts. 



dz £PPP [z] ... 



dz 



Cr 2j7rc^ 



,3 F[Z\[V) r , iQ 



F[z](r) ^ z'Y^ [A'"'''°M/;^(r)I^,^,(r) + ^^^''°W^^(r).W:,(r)]«>2 -Q (^2 ^ „2 ^ 
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Im[z] n+i|z3l 
+ia3 >!: 

6 +i Izjl 
+ia, i 

^ I 

g+i lz,l 

A-ilZil Re[z] 

-i Izjl 

1 -i IzJ 
I ^ 

I 

FIG. 4: Branch cuts of [p^^(r)]°. The branch cuts join the 
integrable poles of [p"'^^(r)]°' at z — ±i\un/v^\ (squares) and 
its zeros at z = iiop, (crosses). 

that < < y < jz^j, as well as on the entire 

interval [— zij+^i], where zi denotes the closest pole to 
the origin. The corresponding branch cuts are charac- 
terized in Fig. m by solid lines joining the zeros of /9^^(r) 
at z = ±ia^ (crosses) and its next integrable pole at 
z = ±z^ (square). Whereas the poles of p^^(r) are in- 
dependent of the position vector r, the points z = ztia^ 
at which it changes sign in between two poles do depend 
on r. 

2. Calculation of E'^Ippp"] 

Knowing the non-analytical structure of the integrand 
F[z](r) [p-'^^(r)] /z^'^-^, the integration contour to be 
used in Eq. ([53]) can be specified. Just as before, the 
circle Cr needs to be deformed in order to apply the 
Cauchy theorem on contours encircling regions where the 
function is entirely analytical. In particular, one cannot 
go through branch cuts as one must remain on the same 
Riemann sheet. An acceptable decomposition under the 
form Cji = [Cg + Cd]{e 0), where each scmi-circlc 
Cg/Cd is further closed by a vertical segment along the 
imaginary axis interrupted by a semi-circle around the 
origin, is displayed in Fig. [5l Note that, as opposed to 
Fig. El no special care needs to be taken around the poles 
at z = ±Zp as they are now integrable {'^ 1/z" with 
< a < 1). The crucial point, however, is that the por- 
tions along the branch cuts will not cancel out as we sum 
the two vertical segments because the integrand (in fact 
[p^^(r)] ) is discontinuous across the branch cuts. 

One may wonder what happens when, as in Fig [SI 
the radius R is such that the original contour Cr goes 
through a branch cut. In fact, the contour Cr defined 



Im[z] h +Z3 CR=Cg+C, 




FIG. 5: Specification of the integration contour for an EDF 
containing fractional powers of the densities. 

through [Cg + Cd] (e ^ 0) in Fig. (i) is well defined 
when a branch cut lies in between Cg and Cd because 
the limit e ^ docs not pose any problem once the value 
of the function on both sides of the cut has been prop- 
erly worked out, (ii) is the contour which has been used 
in actual calculations [J, i, [13 and (iii) might however 
need to be discretized on a rather dense mesh to provide 
converged calculations. 

Note that the deformation of the contour discussed 
above was advocated in Ref. [l| as a remedy to the pathol- 
ogy brought about by branch cuts. In fact, it is rather 
a necessary step to simply define the integration over 
the original circle and obtain the result it provides. As 
detailed below, proceeding to such a deformation of the 
contour does not remove the intrinsic pathological nature 
of MR calculations performed using an EDF containing 
non-integer powers of the density matrices. 

We are now ready to apply the Cauchy theorem along 
the two closed contours appearing in Fig. [5] and then let 
e goes to zero. It is clear that the contributions from the 
vertical portions in between the branch cuts cancel out 
as we add the results from the two closed contours. On 
the other hand, contributions from segments along the 
branch cuts will not cancel out because of the disconti- 
nuity of the integrand across them. 

We consider for illustration (sec Fig. (5]) the situation 
where the contour Cr "hits" the [n -f 1)*'' branch cut 
at z — ±iR; i.e. ctn+i < R < l^n+il- This means that 
the n*'' branch cut is entirely located inside Cr whereas 
the (n -I- 1)*'* one is partially outside the circle of integra- 
tion. For simplicity, and because it is irrelevant to the 
present discussion, we do not calculate the contribution 
£^[ppp"]([— zi, +zi]) from the closest branch cut to the 
origin. Indeed, this one is trickier than the other branch 
cuts because the physical pole at z = lies on that branch 
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Integrable pole 



Branch cut 



cut. All that matters for the present discussion is that 
the branch cut [— 2;i,+Zi] provides a finite contribution 
to the projected energy. In the end, one obtains 



FIG. 6: Zoom on the integration contour Cr obtained as the 
limit of the sum of two disconnected semi-circles. For illus- 
tration, we display a situation where the chosen integration 
contour Cr "hits" the (n -I- l)"* branch cut at z = ±ii?, that 
is, has a radius R such that Qn+i < R< \zn+i\- 



f^[ppp"](i?)-£^[ppp"](hzi,+zi]) = (-1)^ - sin(a7r) 

TT 

which is real and where, for y real. 



dy 



3 F[zy](r) 

niN+l 



p^^y{v) (55) 



1 2 

W7, V 



(56) 
(57) 



C>0 



The above analytical results are explicit enough that wc 
can draw several important conclusions from them. First, 
Eq. ([55)) demonstrates that the PNR energy depends 
on the radius K of the integration contour through the 
boundary of the integral; i.e. the PNR energy is not shift 
invariant. As Cb. goes through a branch cut, the con- 
tribution of that branch cut changes progressively and 
leaves a smoothed step in the PNR energy; see Fig. [T] 
This relates to an unphysical breaking of shift invari- 
ance. Second, there is no discontinuity or divergence as 
Cfl passes through the branch points since the function 
|p^'^(r)|" is integrable at y = jz^j, for all ^. 

The two previous conclusions are at variance with what 
happens for (most of the) EDFs containing only integer 
powers of the densities as recalled in Sect. IIV Bl Indeed, 
a pole crossing the integration provides in this case PNR 
energies with (i) an abrupt step (ii) a divergence if the 
pole is of even order Q. Also, it is important to under- 



line the role played by the regularization of the bilinear 
part of E''^<' put forward in Sect. IIV C 21 If one were to 
use the uncorrected term £ppp , the PNR energy would 
diverge as Cr passes through the branch points. Indeed, 
the integrand in Eq. (j55p would then contain terms over- 
all proportional to {y'^ — |z^p)~^|p^*^(r)|" which is not 
integrable at y = jz^j. 

In any case, the absence of divergence for the regular- 
ized S'^Ippp"] is critical since the associated integrability 
of the pole was used in Ref. Q to assess the meaning- 
fulness of PNR calculations performed with the Gogny 
force. However, and although divergences do constitute 
a dramatic pathology of ill-defined PNR calculations, the 
most profound problem relates rather to the breaking 
of shift invariance of the PNR energy as one changes 
the integration contour. Indeed, the associated spuri- 
ous branch cuts modify the topology of potential energy 
curves as one deforms the system with respect to a col- 
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de'^[ppp"] 
dR 






IzJ R 








R 



No branch cut 



lective degree of freedom. As discussed above, such a 
problem persists for a regularized non-integer power or, 
equivalently, for an effective two-body vertex depending 
on a fractional power of the density. Still, the absence of 
divergence explains why the spurious nature of fractional 
powers of the densities that we focus on here has been 
overlooked so far even more than the pathologies brought 
about by integer powers. 



p.* branch cut 



FIG. 7: Schematic effect of a shift transformation on the PNR 
energy. Top: projected energy E^lppp"'] as a function of R. 
Bottom: same for the derivative of S^lppp""] with respect to 
R. 



In the end, divergences are presently replaced by an- 
other pathological behavior of the PNR energy. To 
isolate such a pattern, let us take the derivative of 
£'^[ppp"]{R) in Eq. ([55]) with respect to the radius of 
integration R. One obtains, for fj. > 1 



dR 



sin(a7r) / d^r F[iR]{r) Ip^'^ir) 



if R e [|z^_i|,q;^] , 
if i? G [a^, Iz^l]. 



(58) 



Because of the non-analytic behavior of |/9^'^(r)|" at 
each branch point, the derivative diverges in Eq. (|58p for 
R ~ \z^\, fi 1. As a result, the projected energy dis- 
plays a kink (non-derivable behavior) as the integration 
circle goes through a branch point or as a branch point 
goes through the integration circle when the system is 
deformed along a collective path. This fact alone is un- 
acceptable for a well-defined projected theory. The corre- 
sponding pattern is schematically displayed in Fig. [7] and 
is observed in realistic calculations as will be discussed 
in Sect. IVl 



and q = p for neutrons and protons, respectively. The 
problematic terms entering the toy Skyrme functional 
(Eqs. [^011^ now take the form 



(59) 



q,q' =p,'. 



+ B^^''°s,(r)-s,,(r) p^{r)l60) 



E. Isospin degree of freedom 

The isospin degree of freedom does not modify any con- 
clusion of the present paper but only complexifies certain 
aspects of the discussion. Still, to provide an idea of the 
modifications brought about by the consideration of both 
protons and neutrons, we now proceed to a restricted set 
of remarks. 

Considering the isospin degree of freedom, one must 
account for the fact that densities, e.g. Pq{v), and single- 
particle wave-functions Lp^[rq) are now labeled with the 
isospin projection quantum number q, where q = n 



S-^P^ ^ /dV^ A'^'^"^ |p,(r)|Vj(r), (61) 

q—p^n 

where the coupling constants A/B characterize terms in 
which the two linear densities involved refer to identi- 
cal/different isospins. Note that neutron-proton pairing 
is not considered. Also, pn{r) is the isoscalar part of the 
matter density. As single-particle states have a definite 
isospin projection, pQ{r) = p„(r) + Pp{r). 

In the present case, both neutron and proton particle 
numbers are restored. Doing so requires to consider two 
gauge angles ifn and ipp for neutrons and protons, respec- 
tively. As a result, PNR energies are obtained through 
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a double integration over the complex plane where the 
corresponding variables are denoted as Zn and Zp. 

As far as the rcgularization of the bilinear part of the 
toy functional, see Sect. IIV C 2| it still leads to the con- 
dition A^^ = —APP and thus only constrains the like- 
particle interaction. Then, one notes that the pseudo 
matrix elements introduced in Eq. (j52p to deal with the 
part of the EDF containing non-integer powers of the 
density matrices now depend on both the neutron z„ 
and proton Zp gauge variables because of the dependence 
on the isoscalar part of the transition local density in 
Eqs. (|59ll6ip . With the pseudo two-body matrix elements 
yPPJl^A'^m Zp] at hand, one can apply the correction for- 
mula of Eq. (43) of Paper II ensuring that the matrix 
elements are now placed underneath the integrals over 
the two gauge angles. 

Once the part of the energy kernel £[zn,Zp] that de- 
pends only on integer powers of the density matrix has 
been regularized, one is left with the spuriosities brought 
by the fractional power of the isoscalar transition den- 
sity [/3g^'(r) -I- pg^''(r)]". The branch cuts of the latter 
are not the same as those seen when dealing with a sin- 
gle particle species. This modifies the analysis but does 
not change the fact that the theory is not satisfactory, 
irrespective of the fine tuning done to define the integra- 
tion contour. As a result, PNR energies cannot be made 
shift invariant and display smooth spurious steps as one 
changes the proton and/or neutron radii of integration 
or deforms the system along a certain degree of freedom. 

V. APPLICATIONS 

We wish to illustrate the analytical results obtained in 
the previous Sections through results of realistic calcu- 
lations. We perform PNR calculations after variation of 
^^O. We use the SLy4 parametrization [3 of the Skyrme 
EDF together with a pairing functional derived from a 
Delta Interaction (DI). The Coulomb exchange part of 
the functional, usually calculated in the Slater approxi- 
mation, is omitted as done in Paper II. The SLy4 Skyrme 
parametrization includes a term of the type S.ppp^''^ which 
is perfectly suited to the present discussion. 

A. Uncorrected calculations 

As explained in Sect. IIII Al traditional PNR calcula- 
tions have been performed using non-diagonal kernels de- 
fined through the prescription £[0, ip\ = £[p^f , k"'^, 
where k, k*] is the single-reference EDF. Figure [8] 
shows the PNR energy £^ obtained in this way for ^^O 
and displayed as a function of quadrupolc deformation. 
The calculation is repeated twice, using 5 and 199 points 
in the discretization of the integrals over the two gauge 
angles. 

One observes that the deformation energy surface ob- 
tained with 5 integration points is smooth and looks 
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FIG. 8: (Color online) Spectrum of poles = for 
protons (top panel) and neutrons (middle panel) as a func- 
tion of quadrupole deformation, which for levels in the vicin- 
ity of the Fermi energy resembles a stretched and slightly 
distorted Nilsson diagram. The dashed red line ai Zq — 1 de- 
notes the radius of the standard integration-contour Rq — 1. 
The bottom panel shows the PNR energy £^ for two differ- 
ent numbers of discretization points in the computation of the 
integrals over the gauge neutron ip„ and proton tpp angles. 



physically reasonable. However, as one increases the 
number of integration points, divergences develop, pre- 
cisely at deformations where a neutron or a proton single- 
particle state crosses the Fermi energy in the underlying 
SR states, i.e. when the associated non-intcgrable branch 
point crosses the unit circle in the complex plane. This 
is consistent with the discussion given in Sect. IIV D 21 for 
the uncorrected SLy4 parametrization. Such divergences 
are at variance with the results obtained in Paper II with 
the SIII parametrization. Indeed, SIII is of specific func- 
tional form such that all the poles at z = are simple 
poles. This is notably due to the fact that the trilincar 
terms entering SIII do not display products of three den- 
sity matrices referring to the same isospin. As explain 
in Sect. IIVB 3[ this property leads to a finite Cauchy 
principle value as the poles cross the integration circle. 
Still, the finite step left in the PNR energy as a 
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FIG. 9: Energy gain from particle number restoration as a 
function of quadrupole deformation for two different numbers 
of discretization points in the computation of the integrals 
over the gauge angles. 



pole/branch cut enters or leaves the integration circle is 
a pathology shared by the calculations performed with 
SLy4 and SIII. Those steps are better visible in Fig. [5] 
which displays the gain from particle number restoration 
with respect to the SR energy (rather than the absolute 
PNR binding energy) using SLy4. Note in passing that 
the reason why the structure around (32 = 0.7 does not 
display a typical step can be understood from the fact 
that two pairs of levels cross the Fermi energy at that 
deformation, as discussed in Paper II. 

By looking carefully, one can observe an interesting 
difference between the steps produced by SIII (see Pa- 
per II) and those obtained presently using SLy4. The 
steps generated by SLy4 are significantly less steep than 
those produced by SIII. This is because, whereas a sharp 
step is generated by an isolated pole leaving or entering 
the integration circle in the case of SIII, which occurs 
over an infinitesimal interval of deformation, it is gener- 
ated by a branch cut leaving or entering the integration 
circle in the case of SLy4, which happens over a finite 
interval of deformation. 



B. Correcting the bilinear part 

The specificity of SLy4 is to contain a term of the type 
£PPP . As discussed in Sect. UVCl one could have hoped 
that regularizing the quadratic part of this term through 
the correction method proposed in Paper I would lead to 
a well-behaved PNR energy; i.e. that the remaining frac- 
tional power of the density would not create any pathol- 
ogy, in particular in view of the fact that the branch point 
becomes integrablc in this case. Of course, it is impor- 
tant to remember that the correction method proposed 
in Paper I relies on solid basis only for terms of the form 



FIG. 10: (Color online) Particle number restored energy £^ 
as a function of quadrupole deformation without and with 
regularization of all bilinear terms in the EDF, including the 
quadratic part of £'''''' ' . Results are shown for two different 
numbers of discretization points in the computation of the 
integrals over the gauge angles. 

EP , with n integer. Thus, regularizing the quadratic 
part of £PPP ^ in this way is purely empirical. 

As a matter of fact, the results displayed in Fig. [10] 
demonstrate that proceeding to such a correction does 
not lead to a well-behaved PNR energy. The integrabil- 
ity of the branch points remaining after regularizing the 
quadratic part of E^pp^'^ is such that all the divergences 
have disappeared. This is a necessary but not sufficient 
condition to obtain a well-behaved PNR energy. Indeed, 
Fig. [n] clearly demonstrates that the spurious steps are 
still present and have in fact not been reduced by regu- 
larizing the bilinear part of Eppp ' . In addition, one ob- 
serves that the corrected results still depend strongly on 
the discretization of the integrals over the gauge angles. 
More precisely, all terms of the energy functional that are 
strictly bilinear have become independent on the number 
of discretization points whereas the term with the extra 
fractional power is not. Considering the experience we 
have gathered about well-behaved PNR energies, such a 
dependence is a fingerprint of a ill-defined PNR theory. 

As discussed in Sect. lIVD2l Figs. [TOl and [TT] also show 
that regularizing the quadratic part of Eppp^'^ leads to 
the replacement of divergences by non-derivable points 
in the PNR potential energy curve. Indeed, kinks are 
clearly visible at the deformation where the divergences 
appeared before applying the correction method. Using 
more mesh points for Qiq, and ipn, one could resolve 
even better the non-derivable character of the energy as a 
branch point passes through the integration circle. This 
pattern relates directly to the analytical result obtained 
in Eq. 

Finally, note that it is a particularity of the SLy4 
interaction complemented with the pairing interaction 
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FIG. 11: (Color online) Energy gain from PNR as a function 
of quadrupole deformation without and with regularization 
of all bilinear terms in the EDF, including the quadratic part 
of S""" . Resuhs are shown for two different numbers of 
discretization points in the computation of the integrals over 
the gauge angles. 



chosen here that the combined correction of all density- 
independent terms is always very small in ^^O, often even 
difficult to resolve on the plots. 



C. Shift transformation 

The finite steps that arise in the deformation energy 
surface are a reminiscence of the violation of the shift 
invariance of the PNR energy. Such a violation is unam- 
biguously demonstrated by varying the radius of the in- 
tegration contour in Eq. ([55]) : i.e. by computing Eq. ([55]) 
as a function of R. 

The upper panel of Fig. [12] shows the PNR energy of 
^^O at a deformation Q20 = 600 fm^, obtained using 
the SLy4 parametrization. The energy is displayed as 
a function of the radius of the integration contour used 
to restore the proton number. The radius for the neu- 
trons is Rn = 1 in all cases. The calculation is performed 
with and without a regularization of the bilinear part of 
the functional and for two different numbers of integra- 
tion points (taken to be the same for protons and neu- 
trons). Finally, the bottom panel of Fig. [T^ shows the 
same quantity obtained from the SIII parametrization at 
a quadrupole deformation Q2Q — 500 fm^. 

The upper panel of Fig. [T^] confirms that, even after 
regularizing the bilinear part of £ppp the PNR en- 
ergy is not invariant under shift transformation. Even 
though the correction method does remove the diver- 
gence, it does not eliminate the shaped steps as the in- 
tegration contour goes through a branch cut. In addi- 
tion, both the corrected and uncorrected PNR energies 
depends strongly on the discretization of the integrals. 
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FIG. 12: (Color online) Particle-number restored energy £^ 
as a function of the radius Rp of the contour chosen to re- 
store proton number {Rn = 1) and for two different num- 
bers of discretization points in the computation of the inte- 
grals over the gauge angles. Results are shown before and 
after regularization of the bilinear part of the EDF. Upper 
panel: at a prolate quadrupole deformation Q20 = 600 fm^ 
using the SLy4 parametrization. Bottom panel: at a pro- 
late quadrupole deformation Q20 = 500 fm^ using the SIII 
parametrization. The corrected SIII curve is independent on 
the number of discretization point; hence, only one curve is 
shown. The left scale shows the absolute value of the binding 
energy whereas the right scale shows the energy gain from 
symmetry restoration. 



Again, those two features are entirely due to the term in 
the functional depending on a non-integer power of the 
density. After regularization, all terms that are strictly 
bilinear become shift invariant. For comparison, the bot- 
tom panel of Fig. [12] shows the PNR energy obtained 
with SIII in Paper II. We recall that SIII contains only 
linear, bilinear and trilincar terms which are such that 
all poles at z = iz^ are of order one. The corresponding 
PNR energy is, after regularization, independent on the 
contour and the number of discretization points with a 
numerical precision better than 1 keV. When restoring 
the particle number that the SR-EDF calculation was 
constrained to, the finite spurious contributions are the 
smallest when using a circle radius close to i? = 1 for the 
reasons outlined in Paper II. Consequently, the corrected 
value is rather close to the uncorrected one in such a case. 
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It is fortuitous that for the deformation Q20 ~ 500 fm^ 
in "'^^O and when using SLy4 the combined correction of 
ah density-independent terms is very small, such that 
corrected and uncorrected curves are close at very small 
values of Rp in Fig. 1121 and even cannot be distinguished 
within the resolution of the plot for larger Rp shown. 

Just as for the deformation energy curve as a func- 
tion of quadrupole deformation, one observes, by com- 
paring the two panels of Fig. [121 that the steps generated 
by SLy4 are significantly less steep than those produced 
by SIII before correction (calculated in both cases with 
enough integration points to resolve them). This is due 
to the fact that the steps are generated by a single pole 
leaving or entering the integration circle in the case of 
SIII, which occurs over an infinitesimal variation of i?p, 
whereas they are generated by a branch cut leaving or en- 
tering the integration contour in the case of SLy4, which 
happens over a finite interval of variation of Rp . 

Just as for the behavior of the deformation energy 
curve as a function of quadrupole deformation, the curves 
obtained with 199 integration points in the upper panel 
of Fig. [T2I show that the divergences seen before regular- 
izing the quadratic part of E^^p have been replaced by 
cusps. Using more mesh points for Rp and ipp, one could 
resolve even better the non-derivable character of the 
PNR energy as the integration circle passes the branch 
points. This is a direct illustration of the analytical re- 
sult obtained in Eq. fSS]) and is schematically displayed 
in Fig. El 

An important byproduct of the previous result is that 
they invalidate PNR calculations performed using a fully 
antisymmetrized two-body interaction that depends on 
the medium through a fractional power of the density, e.g. 
the Gogny interaction. The problem was further circum- 
vented in Ref. Q by using the projected density in place 
of the transition density in the density-dependent term 
of the Gogny interaction. However, such a procedure sin- 
gles out one density factor in the energy kernel in a way 
that seems highly arbitrary and not easily extendable to 
more involved EDFs. In addition, such a prescription of 
using the correlated density into the density-dependent 
term of the effective vertex leads to unsatisfactory results 
for other multi-reference calculations; e.g. calculations in- 
cluding parity restoration and configuration mixing along 
the octupolc degree of freedom [l^ . 



VI. SUMMARY AND CONCLUSIONS 

In Ref. ill , pathologies of calculations aiming at restor- 
ing particle number and performed within the Energy 
Density Functional (EDF) framework have been high- 
lighted. In Ref. 0, the first paper of the present 
series, we demonstrated that such pathologies are in 
fact shared by all multi-reference (MR) calculations, 
i.e. symmetry restoration and/or Generator Coordinate 
Method (GCM)-based configuration mixing calculations, 
performed within the EDF framework. In Ref. Q, a for- 



mal and practical solution that applies (i) to any symme- 
try restoration and/or GCM-based configuration mixing 
calculation (ii) to EDFs depending only on integer pow- 
ers of the density matrices, was proposed. In Ref. 
the second paper of the present series, the regulariza- 
tion method was applied to Particle Number Restoration 
(PNR) calculations using an energy functional that de- 
pends only on integer powers of the density matrices; e.g. 
which contains linear, bilinear and trilinear terms. 

The limitation of the correction method proposed in 
Ref. [5| to energy functionals depending on integer powers 
of the density matrices is a critical feature as most func- 
tionals found in the literature contain non-integer powers 
of the (normal) density matrix, both in the functional 
modeling the strong interaction and in the functional 
modeling the Coulomb interaction, due to the Slater ap- 
proximation to the exchange term Q. Such non-integer 
powers of the density matrices pose difficulties which 
go beyond those posed by integer powers: as transition 
densities are complex, taking their non-integer powers 
amounts to dealing with a multivalued function on the 
complex plane. This makes the analysis of the associated 
pathologies more involved. 

In the present paper, the third of the series, the viabil- 
ity of non-integer powers of the density matrices has been 
addressed, building upon the analysis already carried out 
in Ref. [l|. First, we proposed to reduce the pathologi- 
cal character of terms depending on a non-integer power 
of the density matrices by regularizing the fraction that 
relates to the integer part of the exponent, using the 
method proposed in Ref. [Hi- This amounts to scaling 
down the extent of the problem to the one potentially 
encountered using a fully antisymmetrized effective in- 
teraction depending further on a fractional power of the 
density; e.g. the Gogny force. Second, we discussed in 
detail the spurious character of the remaining fractional 
power of the density (matrix). Both through analytical 
derivations and numerical applications (using the SLy4 
Skyrme parametrization), we demonstrated that regu- 
larizing the fraction related to the integer part of the 
exponent does remove divergences in the particle num- 
ber restored energy but replace them by cusps which are 
as unphysical as the original divergences. In addition, 
the spurious steps in the PNR energy and the related 
breaking of shift invariancc prevail. Such results thus 
invalidate PNR calculations performed using a fully an- 
tisymmetrized two-body interaction that depends on the 
medium through a fractional power of the density. 

Eventually, and because we do not see any well-defined 
basis to correct the corresponding pathologies, we con- 
clude at this point that non-integer powers of the den- 
sity matrices are not viable and should be avoided in 
the first place when constructing nuclear energy density 
functionals to be used in MR-EDF calculations in the fu- 
ture. However, one will have to restrict the form to rather 
low integer orders in the density matrices. For example, 
the EDF recently proposed by Baldo et al. [l^l includes 
terms up to fifth power in the total density p(r), which 
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lead to self-interaction terms [2l| that will require a regu- 
larization containing quadruple sums over single-particle 
states, which will be too costly in realistic calculations. 

Let us make an additional comment regarding the dras- 
tic conclusion to discard non-integer powers of the den- 
sity matrices altogether. On the one hand, integer powers 
of the density matrices appear naturally when construct- 
ing the EDF through ab-initio calculations, e.g. through 
many-body perturbation theory. On the other hand, 
non-integer powers of the density matrices, if not intro- 
duced merely on phenomenological grounds, do often, if 
not always, result from interpreting integrals over mo- 
menta up to kp providing the infinite matter equation 
of state with contributions of the kind kp as density- 
dependent term through the use of kp ~ pl/^ rp^. 

ans- 

ported to finite nuclei, where the latter relationship has 
no rigorous basis, through some version of the local den- 
sity approximation, this leads to an EDF that contains 
terms of the form p"/^. Although such a constructive 
procedure of the nuclear EDF does not lead to particular 
problems in single reference (SR) calculations, it does so 
when this procedure is extended to MR calculations as 



even the local part of the scalar-isoscalar transition den- 
sity matrix is complex, stretching one step too far the 
above procedure proceeding through infinite matter and 
the use of kp <^ p^/'^. Finally, there are both practi- 
cal reasons and formal motivations to conclude that (i) 
non-integer powers of the density (matrix) are not viable 
in (multi-reference) EDF calculations (ii) paramcteriza- 
tions making only use of integer powers of the densities 
need to be constructed in the very near future. Last but 
not least, note that such a conclusion actually extends to 
any form of the EDF that generates branch cuts when 
continued over the complex. 
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